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There is great interest in numerical relativity simulations involving matter due to the likelihood 
that binary compact objects involving neutron stars will be detected by gravitational wave obser- 
vatories in the coming years, as well as to the possibility that binary compact object mergers could 
explain short-duration gamma-ray bursts. We present a code designed for simulations of hydrody- 
namics coupled to the Einstein field equations targeted toward such applications. This code has 
recently been used to study eccentric mergers of black hole-neutron star binaries. We evolve the fluid 
conservatively using high-resolution shock-capturing methods, while the field equations are solved 
in the generalized-harmonic formulation with finite differences. In order to resolve the various scales 
that may arise, we use adaptive mesh refinement (AMR) with grid hierarchies based on truncation 
error estimates. A noteworthy feature of this code is the implementation of the flux correction al- 
gorithm of Berger and Colella to ensure that the conservative nature of fluid advection is respected 
across AMR boundaries. We present various tests to compare the performance of different limiters 
and flux calculation methods, as well as to demonstrate the utility of AMR flux corrections. 



I. INTRODUCTION 

The interface between strong field gravity and mat- 
ter dynamics promises to be one of the important fron- 
tiers in the coming years. A new generation of gravita- 
tional wave detectors (LIGO Q, GEO 0,TAMA||, and 
VIRGO [||]) are now operational, and within the next few 
years are expected to reach sensitivities that will allow 
observations of the Universe in gravitational radiation for 
the first time. The prime targets of these observations 
are compact object (CO) binaries composed of combi- 
nations of black holes (BHs) and neutron stars (NSs). 
Modeling of such sources is a crucial ingredient to real- 
ize the promise of gravitational wave astronomy. Even if 
an event is detected with a high signal-to-noise ratio, re- 
constructing the dynamics of the system that produced 
the signal cannot be done directly but instead will re- 
quire template banks of theoretical waveforms informed 
by numerical simulations. 

Compact object mergers involving NSs are expected 
to be significant sources of not only gravitational radia- 
tion, but also possible progenitors for short-gamma-ray 
bursts (SGRBs) 043 and other electromagnetic and neu- 
trino counterparts [8|. Efforts are already underway to 
use potential gravitational wave sources as triggers for 
searches for electromagnetic transients @, [To[ • Observa- 
tions would help constrain evolutionary models for the 
progenitor stars and their environments. Perhaps most 
intriguingly, the observations would give clues to the 
equation of state (EOS) of matter at nuclear densities (as 
in NS interiors), which cannot be probed in laboratories 
on Earth and is not fully understood at the theoretical 
level (for a broad discussion see for example [HI). The 
reason that the gravitational wave signature could con- 
tain information about the matter EOS (and other details 
about the internal structure of neutron stars) is that the 
EOS in general has a significant effect on the bulk motion 
of matter, and it is this bulk motion that is the mecha- 
nism by which gravitational waves are produced. Several 



studies to date have looked into this issue, suggesting the 
imprint of the EOS on the gravitational waves will be 
strong enough to detect |12h23| (though, in some cases, 
the expected frequencies are higher than the range to 
which the current generation of ground-based detectors 
are most sensitive, thus limiting the information which 
can be extracted). While CO binaries containing NSs 
are a particularly interesting class of sources involving 
general relativistic (GR) hydrodynamics, they are by no 
means the only such systems. Examples of additional 
systems that have already been considered include BH 
accretion tori .24-[27| and NS-white dwarf mergers pjj . 

Thoroughly modeling systems like those described 
above would require evolution of the spacetime, the pho- 
ton and neutrino radiation fields, and the magnetized, 
relativistic fluid. Even a minimalistic treatment, with 
the Einstein equations coupled to the equations of rela- 
tivistic hydrodynamics, represents a complex, nonlinear 
system of partial differential equations. Numerical simu- 
lations are thus essential for exploring such strong field, 
dynamical systems. There is a long history of adapting 
successful techniques for simulating Newtonian hydrody- 
namics to relativistic and general relativistic fluids which 
we will not attempt to summarize (see [2^ for a review 
of general relativistic hydrodynamics). Instead, we will 
briefly attempt to place the code described in the present 
paper in the context of other recent codes developed for 
fluids on evolving spacetimes. R 

Several of these codes 1334301 solve the field equations 
in the BSSN formulation [37|,438| . The remainder HUES 
use the generalized-harmonic formulation [4l|, [42[ which 



1 Note that our focus is restricted to codes which handle dynam- 
ically evolving gravitational fields. Such codes, however, fre- 
quently owe much to earlier, fixed-background evolution codes 
(see [23|). In addition, advancements such as GR- hydro with 
multipatch grids l3ll and with GPUs l3ll have recently been 
made with fixed-background codes. 
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we also employ; unlike our code, however, these groups 
convert to a fully first-order formulation (43|. Most 
groups use finite-difference methods for the metric evolu- 
tion and a conservative, high-resolution shock-capturing 
(HRSC) scheme for the hydro evolution; these unigrid 
algorithms are then interfaced with some sort of adap- 
tive mesh refinement (AMR). A notable exception for 
the metric evolution is (39j . which employs pseudospec- 
tral methods for the metric and then interpolates to a 
finite- volume grid for the fluid. 

Some groups have implemented the MHD equations 
in full GR; since these codes all make use of conserva- 
tive HRSC methods, they may be principally differenti- 
ated by how they meet the challenge of preserving the 
V • B = constraint. (A straightforward finite-difference 
evolution of the magnetic field would generically lead to 
magnetic monopoles and, hence, unphysical behavior.) 
WhiskyMHD employs constrained transport [34[ for this 
purpose, which preserves the constraint to machine ac- 
curacy, whereas the code of 44[ uses hyperbolic diver- 
gence cleaning. Constrained transport, however, requires 
special interpolation at refinement level boundaries in or- 
der to preserve the constraint. The Illinois group found 
that a vector-potential formulation of the MHD equa- 
tions works well when coupled to AMR (4f|. This is be- 
cause the constraint is preserved by construction with the 
vector-potential, even with the restriction and prolonga- 
tion operations of AMR (see also [46| for a thorough ex- 
amination of the electromagnetic gauge condition) . Stud- 
ies indicate that magnetic fields do not significantly affect 
the gravitational dynamics of CO mergers (see e.g. 44j). 



but they could be critical for understanding EM counter- 
parts including the possible formation of a SGRB engine. 
A new method to treat the MHD equations was recently 
presented in 47[, where ideal MHD is used in high mat- 
ter density regions {e.g. inside a NS), while the force-free 
approximation is used elsewhere {e.g. the magnetosphere 
of a NS). The authors applied the method to study the 
collapse of magnetized hypermassive NSs (which could 
be formed via binary NS mergers) and suggested that 
intense EM outbursts could accompany such events. 

Besides MHD, the other major advances in the phys- 
ical model for numerical relativity codes have been in 
the arena of microphysics. While the r = 2 EOS was 
the community standard for quite some time, most codes 
now allow for a nuclear theory-based EOS [48|,|49| and/or 
use various parametrized, piecewise polytropic EOSs in- 
spired by the range of plausible nuclear EOSs [5(| HH . 
These advances in EOS description primarily affect the 
cold NS structure, but the group developing the SACRA 
code has also begun to account for neutrino transport 
via a simplified leakage scheme fl2l . l52j . The same group 
has also made available a formulation for a more accurate 
truncated moment scheme with a variable Eddington fac- 
tor closure [53j , which shows much promise for numerical 
relativity simulations with neutrino physics beyond the 
leakage approximation. 

Another category of GR hydrodynamics codes employs 



the conformal-flatness approximation, which is particu- 
larly useful when supernova simulations are the target 
application. An example is CoCoNuT/ VERTEX, which 
incorporates relativistic hydrodynamics, conformally flat 
gravity, and ray- by-ray neutrino transport [Hij . The code 
of (55j employs a similar scheme for hydrodynamics and 
gravity but adds a test magnetic field; this code has been 
used to study the magnetorotational instability in super- 
novae. 

Newtonian (and semi- Newtonian) l57j , conformally 
flat [H, H^) , and fixed- background (601 SPH codes rep- 
resent an important, orthogonal approach to studying 
CO interactions. SPH has an advantage over Eulerian 
schemes when a large range of spatial scales is involved. 
Such a situation may arise in CO mergers when material 
is stripped from a star in a tidal interaction and forms 
an extended tail. On the other hand, Eulerian codes are 
the standard approach when strong shocks are present, 
as would arise in binary NS mergers or disk circulariza- 
tion. (Recent progress has been made, however, in ap- 
plying SPH to situations with relativistic shocks [6l|.) In 
addition, SPH has not (to our knowledge) yet been cou- 
pled to a code which evolves the full Einstein equations. 
Nonetheless, comparisons between Eulerian and SPH re- 
sults could prove very useful on a problem-by-problem 
basis to characterize the errors in both methods. 

Though current efforts in GR simulations involving 
matter tend to focus on increasingly complex physical 
models, there remain many unanswered questions in the 
astrophysics of compact objects that can be addressed 
with a code which solves the Einstein equations coupled 
to perfect fluid hydrodynamics. We have thus focused our 
code development on hydrodynamics in full GR, while 
maintaining a flexible infrastructure to accommodate ad- 
ditional physics modules in the future. We evolve the 
field equations in the generalized-harmonic formulation 
using finite differences. The fluid is evolved conserva- 
tively using one of several different shock-capturing tech- 
niques we test here. We have also implemented the hy- 
drodynamical equations in a manner that is independent 
of EOS. We make use of AMR by dynamically adapting 
the mesh refinement hierarchy based on truncation er- 
ror estimates of a select number of the evolved variables. 
We also utilize Berger and Colella [62| style flux correc- 
tions (also known as "rcfluxing") in order to make the 
use of AMR compatible with the conservative nature of 
the hydrodynamic equations. Though AMR flux correc 
tions have been implemented in other astrophysical hy 



drodynamics codes (such as Athena [631 ] ■ CASTRO [64j . 
Enzo (65[, and FLASH (6^|), to our knowledge this algo- 
rithm has not been used previously for hydrodynamics 
simulations in full general relativityo A further notewor- 
thy feature of our implementation is that we store cor- 



2 Note that "flux correction" here refers to the enforcement of con- 
servation at AMR boundaries, not the recalculation of fluxes with 
a more dissipative scheme to preserve stability as in Athena |67|| . 
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rections to the corresponding fluid quantity integrated in 
the volume of a given cell instead of the flux, allowing 
for easy implementation within a computational infras- 
tructure that supports cell-centered but not face-centered 
distributed data structures. The code described here has 
recently been applied to studying BH-NS mergers with 
eccentricity as may arise in dense stellar systems such as 
galactic nuclear clusters and globular clusters [H, Hil • 

In the remainder of this paper we outline our compu- 
tational methodology for simulating hydrodynamics cou- 
pled to the Einstein field equations and describe tests of 
this methodology. In Sec. [TT] we review the generalized- 
harmonic approach to solving the field equations and 
present our methods for conservatively evolving a per- 
fect fluid coupled to gravity, including our method for 
inverting the conserved quantities to obtain the primi- 
tive fluid variables and the implementation of flux correc- 
tions to enforce the conservation of fluid quantities across 
AMR boundaries. In Sec. Mil we present simulation re- 
sults which test these methods, highlight the strengths 
and weaknesses of various shock -capturing techniques, 
and demonstrate the utility of the flux correction algo- 
rithm. 



II. COMPUTATIONAL METHODOLOGY 

In this section we begin by explaining the basic equa- 
tions and variables we use to numerically evolve the Ein- 
stein equations in Sec. Ill Al and then discuss the conserva- 
tive formulation of the hydrodynamics equations that we 
use in Sec. Ill Bl The evolution of conserved fluid variables 
necessitates an algorithm for inverting these quantities to 
obtain the primitive fluid variables which we present in 
Sec. Ill CI Finally in Sec. Ill Dl we present the details of 
our algorithm for AMR with flux corrections. 

A. Solution of the Einstein equations 

We solve the field equations in the generalized- 
harmonic formulation [4lL |42| where we fix the coordi- 
nate degrees of freedom by specifying the evolution of 
the source functions H a := Dx a . In this form the evo- 
lution equation for the metric, g ab , becomes manifestly 
hyperbolic: 

g cd d c d d g ab + d b g cd d c g ad + d a g cd d c g bd 

+ 2 #(a,fc) - 2HdT ab + 2T db v ca 

= -8n(2T ab - g ab T) (1) 

where T^ c is the Christoffcl symbol, T ab is the stress- 
energy tensor, and T is its trace. We evolve the metric, 
the source functions, and their respective time derivatives 
using fourth-order Runge-Kutta where the spatial deriva- 
tives are calculated using fourth-order accurate finite- 
difference techniques. In other words, we have reduced 
the evolution equations to first order in time so that there 



are 28 "fundamental" variables {gab,H a ,d t gab,d t H a }, 
but we directly discretize all first and second spatial gra- 
dients without the introduction of additional auxiliary 
variables. 

Analytically one can show [7(J that if one begins with 
initial data that satisfies the Hamiltonian and momen- 
tum constraints, initially set H a = Ox a , and then evolve 
the metric according to ([1]) and the source functions ac- 
cording to some specified differential equations, then the 
constraint equation H a — Ox a = will be satisfied for 
all time. Numerically this statement will only be true to 
within truncation error, which can grow exponentially in 
black hole space times; to prevent this we add constraint 
damping terms as in [7ll. |72|. In practice, ensuring that 
H a — Dx a is converging to zero for a given numerical sim- 
ulation run at different resolutions provides an excellent 
check that the numerical solution is indeed converging to 
a solution of the field equations. 

As described in [42j , the computational grid we use is 
compactified so as to include spatial infinity. This way 
we can impose boundary conditions on the metric simply 
by requiring that it be Minkowski. However we evolve 
the metric of the uncompactificd coordinates since the 
compactified metric is singular at spatial infinity. 

B. Conservative Hydrodynamics 

Coupled to gravity we consider a perfect fluid with 
stress-energy tensor 

T ab = phu a u b + g ab P , (2) 

where h := 1 + P/ p + e is the specific enthalpy and u a is 
the four- velocity of the fluid element. The intrinsic fluid 
quantities p, the rest-mass density; P, the pressure; and 
e, the specific energy are defined in the comoving frame 
of the fluid element. The equations of hydrod yna mics are 
then written in conservative form as follows [73j: 

d t D + d i (Dv i )=0 (3) 
d t S a + di (y/^fjl**) = \v~gT bc d a g bc (4) 

where v 1 is the coordinate velocity, g is the determinant of 
the metric, and the index i runs over spatial coordinates 
only. Note that (0| explicitly contains the time derivative 
of the metric for index a = t. The conserved variables D 
and S a are defined as follows: 

D := ^f^gpu 1 (5) 
S a := ^f~gT l a (6) 

where D is simply the time component of the matter 4- 
current0. 



3 In some implementations of the GR (magneto)hydrodynamic 
equations, see for e.g. |74| , the analog of St in JBJl that is evolved 
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In some situations wc wish to perform axisymmctric 
simulations where we use the symmetry to reduce the 
computational domain to two dimensions. We do this 
using a modification of the Cartoon method [75[ as de- 
scribed in [43 |. where we take the x-axis as the axis of 
symmetry, and only evolve the z = slice of the space- 
time. For the hydrodynamics this means that effectively 
each fluid cell becomes a cylindrical shell, and we use the 
fact that the Lie derivative of the fluid fields with respect 
to the axisymmetric killing vector are zero to rewrite the 
coordinate divergences in the above equations as 

d^Dv 1 ) = d x {Dv x ) + 2d y 2(yD V y) (7) 

and similarly for di(y/—gT z a ) for the t and x components. 
For the y component there is an additional source term 

diiV^gT'y) = d x (V=gT%) + 2^,,. : //N —/,;'. - 

(Sy+V^P/2/). (8) 

By writing the y flux contribution in terms of d y i we 
ensure that when we discretize our evolution will be con- 
servative with respect to the cylindrical shell volume el- 
ement. We choose a special form for the equation for 
S z : 

d t S z + d x {s/—gT\) + -^tfJ—gTl) = 0, (9) 

since in axisymmetry the quantity yS z is exactly con- 
served (that is, it has no source term). 

The conservative evolution system is solved numeri- 
cally using HRSC schemes. We briefly summarize the dif- 
ferent methods we have implemented and test in this pa- 
per, though the references should be consulted for more 
complete details. For calculating intercell fluxes we have 
implemented HLL [76| , the Roe solver [77j , and the Mar- 
quina flux [78| method. The HLL method is straight- 
forward to implement since it does not require the spec- 
tral decomposition of the flux Jacobian and is based on 
estimates for the largest and smallest signal velocities. 
The Roe solver works by solving the linearized Riemann 
problem obtained using the flux Jacobian at each cell in- 
terface (using the so-called Roe average of the left and 
right states). The Marquina flux method is an extension 
of this idea that avoids the artificial intermediate state 
and switches to a more viscous local Lax-Fricdrich-type 
method from [79| when the characteristic speeds change 



has the rest-mass density D subtracted off. This could provide 
improved results in situations where the rest-mass density is or- 
ders of magnitude larger than the internal or magnetic energy, 
and accuracy in these latter quantities is important. Though 
we have not explored this alternative, in the scenarios studied 
here (in particular since we are not looking at the behavior of 
magnetic fields) the added effect of a small amount of internal 
relative to rest energy on the dynamics of the fluid or metric 
will be negligible, and we expect either definition of St to give 
comparable accuracy results here. 



sign across the interface. Since the latter two methods 
require the spectral decomposition of the flux Jacobian, 
we give it for our particular choice of conserved variables 
in the Appendix. For reconstructing fluid primitive vari- 
ables at cell faces we have implemented MC and min- 
mod 0, PPM M B and WEN05 0, all of which 
may be used interchangeably with any flux method. MC 
and minmod are both slope limiter methods that reduce 
to linear reconstruction for smooth flows. Minmod is 
the more diffusive of the two. In comparison, PPM and 
WEN05 are higher-order reconstruction methods. PPM 
is based on parabolic reconstruction with modifications 
to handle contact discontinuities, avoid spurious oscilla- 
tions from shocks by reducing order, and impose mono- 
tonicity. WEN05 combines three different three-point 
stencils with weights that are determined by a measure 
of the smoothness of the quantity being reconstructed. 
The specific fluid quantities that we reconstruct on the 
cell faces are p, u, and WU l , where u := pe, W is the 
Lorentz factor between the local fluid element and an ob- 
server normal to the constant t hypcrsurfaces, and U l is 
the Eulcrian velocity (the explicit form of which is given 
in the following section). We choose to reconstruct WU l 
instead of simply U l since any finite value of this quantity 
corresponds to a subluminal velocity. 

The fluid is evolved in time using second-order Runge- 
Kutta. Since the fluid is evolved in tandem with the 
metric, the first and second substeps of the fluid Runge- 
Kutta time step are chosen to coincide with the first and 
third substeps of the metric time step. Since the spa- 
tial discretization of the fluid equations that we use is 
only second-order we choose to use second-order time 
stepping for the hydrodynamics and we have not yet 
experimented with higher-order methods. Wc still use 
fourth-order Runge-Kutta for nonvacuum metric evolu- 
tion (even though for evolutions with matter the overall 
convergence rate will be no greater than second-order) 
both for convenience and because in vacuum dominated 
regions we may expect some improvement in accuracy. 
For general rclativistic hydrodynamics wc evolve the fluid 
on a finite subset (though the majority) of the total grid 
(which as mentioned extends to spatial infinity through 
our use of compactified coordinates), and at the outer 
boundary for the fluid we impose an outflow condition. 

Finally, as is common practice for this method of sim- 
ulating hydrodynamics, we require that the fluid density 
never drop below a certain threshold, adding a so-called 
numerical atmosphere. We give this numerical atmo- 
sphere a spatial dependence that makes it less dense ap- 
proaching the boundaries^] and choose a maximum value 



4 In particular we use the reconstruction parameters presented 
in \8% . 

5 Specifically, we perform reconstruction with the stencils and 
weights presented in Section A2 of [84j . 

6 Specifically we let p a tm(x c , y c , z v ) = pcos 2 (x c ) cos 2 (j/ c ) cos 2 (z c ) 
where p is a constant, and (x c ,y c , z c ) are the compactified coor- 
dinates which range from -1 to 1. 
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that makes it dynamically negligible (typically at least 10 
orders of magnitude below the maximum density). The 
atmosphere is initialized using a cold equation of state 
(e.g. a polytropic equation of state). 



C. Primitive Inversion 

The set of hydrodynamical equations is closed by an 
EOS of the form P = P(p, e). While the conserved vari- 
ables S a and D are simply expressed in terms of fluid 
primitive variables (p, P, e, and v l ) and the metric, the 
reverse is not true. This necessitates a numerical inver- 
sion to obtain the primitive variables following each up- 
date of the conserved variables. The method we use is 
similar to the one used in (85j for spherical symmetry. 
First, we decompose the 4-dimensional metric into the 
usual ADM space plus time form 

ds 2 = g a f>dx a dx b 

j l3 (dx i + pdt)(dx j + f3 j dt) (10) 



-a 2 dt 2 



where 7^ is the spatial metric, a the lapse function and 
P 1 the shift vector. Then, from the metric and conserved 
variables we construct two quantities, 

S 2 
E 



fiSiSj = 1 H 2 W 2 (W 2 



-g(HW 2 



1) 

-P) 



(11) 
(12) 



where H := ph and 7 is the determinant of the spatial 
metric. We reduce the problem of calculating the prim- 
itive fluid variables from the metric and conserved vari- 
ables to a one-dimensional root problem, where we begin 
with a guess for H and iteratively converge to the correct 
value such that f(H) = for some function. From ([12")) 
we can choose 



1(H) = Ej 4— g 



HW 2 



P. 



(13) 



Note that given the metric and conserved variables, f(H) 
is only a function of H, and can be computed as follows. 
First, calculate W 2 = (1 + Vl + 4A)/2 where 



A := 



W 2 (W 2 - 1) . 



7# 2 

Then compute p and e from 
and 

e = -H(W 2 - l)/p + WE /(Da) - 1 



(14) 



(15) 



(16) 



respectively. Once p and e are known, P can be obtained 
from the equation of state, and then f(H) above. An 
iterative procedure for solving f(H) = 0, where f(H) 
is calculated as just described, thus gives the primitive 
variables p. P. and e. The three-velocity can then be 
computed from 



U' 



(17) 



where the Eulerian velocity U l is related to the grid 
three-velocity through U l = (v l + j3 l )/a. This inversion 
scheme is implemented so as to allow any EOS of the 
form P = P(p, e); thus, T-law, piecewise polytrope, and 
tabular equations of state such as the finite-temperature 
EOS of Shen et al. [8(| H3] (for a given electron fraction 
Y e ) are all supported. 

In practice we solve for f(H) = numerically using 
Brent's method |88| . which does not require knowledge 
of derivatives and is guaranteed to converge for any con- 
tinuous equation of state as long as one begins with a 
brackelQ around the correct solution. This can be useful 
when dealing with equations of state interpolated from 
tabulated values. One can avoid losing accuracy in the ul- 
trarclativistic and nonrelativistic limit by Taylor expand- 
ing the above inversion formulas (see |85|). for example, 
in 1/A and A, respectively. We have implemented such 
expansions in our primitive inversion algorithm, though 
we have not yet made any significant study of the inver- 
sion calculation in these regimes. 

In some cases the conserved variables will, due to nu- 
merical inaccuracies, evolve to a state that does not cor- 
respond to any physical values for the primitive variables. 
This causes the inversion procedure to fail. This can hap- 
pen in very low density regions that are not dynamically 
important but still must be addressed. We handle such 
situations using a method similar to that of [73j by ig- 
noring the value of St and instead requiring the fluid to 
satisfy a cold equation of state. 



D. AMR with flux corrections 

Many of the problems we arc interested in applying 
this code to involve a range of length scales, and in many 
cases we expect the small length scale features not to 
be volume filling, for example the individual compact 
objects in binary mergers. Such scenarios can be effi- 
ciently resolved with Berger and Oliger style adaptive 
mesh refinement [8{| ■ A description of the variant of the 
algorithm we use can be found in [90j ; here we mention 
some particulars to this implementation, and give a de- 
tailed description of the extension to ensure conservation 
across refinement boundaries. 

The computational domain is decomposed into a hi- 
erarchy of uniform meshes, where finer (child) meshes 
are entirely contained within coarser (parent) meshes. 
The hierarchy is constructed using (primarily) trunca- 
tion error (TE) estimates, which are computed within 



7 The initial bracket for the root finding is chosen by first checking 
if [Ho/(l + 5), Ho(l + &)], where Ho is the value of H computed 
for the primitive variables at the previous time step and 8 > 
is a parameter we take to be 0.4, is a valid bracket around the 
zero of f(H). If it is not, as a failsafe we try successively larger 
brackets with [H Q /{I + 8) n ,H {l + <5) n ] for n > 2. 
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the Bcrgcr and Oliger time subcycling procedure by com- 
paring the solution obtained on adjacent levels of refine- 
ment before the coarser levels are overwritten with the 
solution from the finer level. Typically we only use the 
TE of the metric variables, since fluid variables in general 
develop discontinuities as well as turbulent features that 
do not follow strict convergence. The layout of the AMR 
hierarchy is then periodically adjusted in order to keep 
the TE below some global threshold. In some situations 
we also require that a region where the fluid density is 
above a certain threshold always be covered by a mini- 
mum amount of resolution. This can be used to ensure, 
for example, that the resolution around a NS does not 
temporarily drop below some level even if the TE of the 
metric variables in the neighborhood of the star becomes 
small. 

When setting the values of the metric variables on the 
AMR boundary of a given child level we interpolate from 
the parent level using third-order interpolation in time 
and fourth-order in space. For the cell-centered vari- 
ables, the outer two cells in each spatial direction (for 
a refinement ratio of 2) on a child level are initially set 
using second order interpolation in time and space from 
the parent level. Following the evolution of the child 
level and flux correction applied to the parent level when 
they are in sync as described below, but before the cell- 
centered values on the child level are injected into the 
parent level, the values in the child boundary cells arc 
reset using first-order conservative (spatial) interpolation 
from the parent level (i.e. the value in the child cell is 
set to be the same as that of the parent cell in which the 
child cell is contained). This ensures that the boundary 
cells on the child level are consistent with the correspond- 
ing flux-corrected cells on the parent level but does not 
affect the order of convergence of the scheme since these 
values are not used in the evolution step. During a regrid 
when adding cells to the domain of a refined level we also 
use first-order conservative interpolation from the over- 
lapping parent level to initialize the values of the fluid 
variables at new cells (fourth-order interpolation is used 
for the metric variables). Note that the actual domain 
that is refined is larger than the volume where the TE 
estimate is above threshold by a given buffer in any direc- 
tion. The buffer size and rcgridding interval are chosen so 
that if change in the region of high TE is associated with 
bulk motion of the solution (e.g. the NS moving through 
the domain), this region will never move by more that 
the size of the buffer between regrids. This ensures that 
new cells (for this kind of flow) are always interpolated 
from regions of the parent that are below the maximum 
TE threshold. Thus, though the interpolation operation 
to initialize new cells is first-order, we find the error it 
introduces is negligible (i.e., below the maximum desired 
TE). 

AMR boundaries require special treatment in conser- 
vative hydrodynamics codes however, since the fluxes 
across the boundary of a fine-grid region will not exactly 
match the corresponding flux calculated on the coarse- 



grid due to differing truncation errors. To enforce con- 
servation, we correct the adjacent coarse grid cells using 
the fine-grid fluxes according to the method of Berger and 
Colella 62|. In the remainder of this section we review 
the algorithm and outline our specific implementation. 

We will concentrate on the evolution of J) on a 3- 
dimensional spatial grid, though the remaining conserved 
fluid quantities are treated the same way, and modifica- 
tion to different numbers of spatial dimensions is trivial. 
Equation ((3]) is evolved numerically at a given resolution 
as 

u i.j.k — 



+(^ + iA*-- F «-iaJ/^ 

+ (*iJ,*+V2 " *lj,fc-l/2)/A* 



(18) 



where Dt* - fe is the volume average of D over the k) 
cell at time t = nAt; F? +1 ,^ ■ k is the flux F x = Dv x 
through the (i + 1/2, j, k) cell face; Ax is the x length of 
each cell and so on for the y and z direction. In prac- 
tice the flux values will be calculated with some HRSC 
technique combined with Runge-Kutta, but the specifics 
are not relevant here. Now consider a situation with two 
sequential levels of refinement, L and L + 1, where level 
L+l has a higher resolution with spatial refinement ratio 
of r in each direction, and its domain is a subset of level 
L. (In practice, we always take r = 2.) Here we focus 
the discussion on a left boundary in the x direction, as 
illustrated in Fig. [TJ boundaries along the right face and 
other coordinate directions are treated in a like manner. 

When evolving according to the Bcrgcr-Oliger algo- 
rithm, after each time step of length At is taken on level 
L, r steps of length At/r are taken on level L+l. Then 
the results obtained on L + 1 are injected into level L 
where the levels overlap i.e., the restriction operation is 
performed conservatively by setting the value in the par- 
ent cell to the (coordinate) volume-weighted average of 
the child cells that make up the parent cell. Now on level 
L, the change in D due to flux going through the cell face 
(j,L + 1/2, jLi ^l) on a timestcp will be 



(19) 



On level L + l, the change in D in one fine-level time step 
due to flux passing through one of the r 2 cell faces that 
make up this same interface is 



6D 



L-\-l,j,k,m 



At/r 
Ax/r 



F? L+1+ i/2,j L+ M L+1+ k(tn + mAt/r). (20) 

for j, k, and m G {0, 1, . . . , r — 1}. Now because 
of truncation error, in general the change in the net 
"mass 5 M L := 5D L V L within the coarse-level cell 



8 For the conserved fluid variable D which we focus on for speci- 
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FIG. 1. A visualization of a refinement level boundary and its treatment in the flux correction algorithm. The top shows 
cells in the x direction on refinement level L while the bottom shows equivalent cells for the L + l refinement level (here the 
refinement ratio is 2). Fluxes are symbolized by arrows. On the bottom level the blue cells ("type B" in the discussion in the 
text) and those to the left on level L + l are boundary cells and will have their values set by interpolation from level L following 
an evolution step on level L. Because of truncation error, subsequent evolution on level L + l will give a flux differing from that 
computed on the parent level L. Consequently, when the new fine grid solution is injected back to the parent level (in cells to 
the right of the red/dotted pattern cell), the solution about the boundary will no longer be consistent with the flux previously 
computed there. The correct this, the fluid quantity in the red/dotted pattern cell is adjusted to exactly compensate for the 
difference in flux computed between the coarse and fine levels. 




a t («L,jL,fci) computed with the coarse-level fluxes 
will not equal the corresponding quantity SMl + i := 
J2j,k, m SD L+i,]M,mV L +i,j,k,m computed with the fine- 
level fluxes, where Vl is the coordinate volume of the 
cell jL, and VL+ij,fc,m is the coordinate volume 
of the cell (il+i, Jl+i +j, fci+i +k). Thus, after the val- 
ues of D on level L + l are injected into level L (in cells 
(«l + 1, 3L-, kh) and to the right in this example), the solu- 
tion on level L will suffer a violation of mass conservation 
proportional to SMl — SMl+i- To restore the conserva- 
tive nature of the algorithm, the idea, described in detail 
below, is to adjust the conservative variable D in the 
cell (jl,Jl,^l) post-injection by an amount to exactly 
compensate for this truncation error induced difference. 

The scheme originally proposed in [62| is to define 
an array that keeps track of a correction to the fluxes 
through cell faces on level L that make up the boundary 
of the evolved cells on level L + l. Consider the case 
where (%l + 1/2, Jl, &l) is such a face. This face-centered 



ficity, the value of the quantity integrated within the volume of 
a cell in fact represents the rest mass in that cell. Through- 
out this section we will therefore use the term 'mass' to refer to 
the value of a conserved fluid variable volume integrated within 
a cell, though for other conserved fluid variables this will not 
correspond to a physical mass. 



flux correction array, SF, is initialized with the inverse 
of the flux in (TH)1) . SF = —F x , , ,„ . ... , , and then 

* " 1>L+±/*,JL+3 ! KL+K ; 

during the course of taking the r time steps on level L + l 
receives corrections from the terms in (|20j) 

SF ^ 5F +^> E K + ,+i/2, jL+ML+ ,+ k it n +rnM/r). 

j,k,m 

(21) 

After the cell values on level L are overwritten by the 
injected values on level L + l where they overlap, the 
cells on level L that abut level L + l though are not 
themselves covered by level L + l cells are corrected with 
the flux stored in SF. 

The way we implement the flux correction algorithm 
is slightly different from this. In particular we wish to 
avoid the added computational complexity of implement- 
ing face-centered grid functions, and therefore we keep 
track of a cell-centered correction. The correction is thus 
also more naturally represented as a correction to the 
fluid quantity integrated within the volume of the cell 
(e.g. for D the rest- mass) rather than a flux. Again 
referring to Fig. [TJ we define the first few cells at the 
boundary of level L + 1 as buffer cells since the calcula- 
tion of flux requires knowledge of the state on both sides 
of the interface. These cells will have their values set 
by interpolation from those in level L. The innermost 
buffer cells for the boundary on level L + 1 we call type 
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B cells (blue cells in the lower half of the figure). These 
are the cells where the level L + 1 contribution to the 
mass correction will be stored. The cell on level L which 
contains the type B cell we will refer to as a type A cell 
(red, dotted-pattern cell). Type A cells are the ones that 
receive mass corrections in this algorithm. For each cell 
on each refinement level we use a bitmask grid function 
that indicates whether the cell is one of the above types 
(A or B), and if so which of the six possible faces (+x, 
—x, +y, —y, +z, —z) abut the boundary. For simplicity 
in the implementation we do not allow grid hierarchies 
where a cell would be both type A and type E0. 

In the following we outline the additional tasks rel- 
ative to the basic Berger-Oliger algorithm that need to 
be performed with our implementation of Berger-Colella. 
Following the spirit of these algorithms, we break down 
the tasks into those the AMR "driver" code implements, 
which do not require knowledge of the specific equations 
being evolved or what physical quantities the variables 
represent, and conversely the "application" steps that 
would need to be implemented by a unigrid application 
code plugging into the driver to become AMR-capablc. 
The driver tasks include the following: 

(i) For the conserved fluid density D, allocate a storage 
grid function to keep track of the associated mass 
correction SM, i.e. the total correction to D within 
the volume of a given cell. 

(ii) Upon initialization set all correction arrays SM to 
zero, and compute the bitmask for the current re- 
finement hierarchy. 

(iii) After any regrid, recompute the bitmask array for 
the new hierarchy. 

(iv) During the stage when buffer cells are set for vari- 
able D at interior boundaries on level L + 1 via 
interpolation from level L, also interpolate the cor- 
rection variable SM , where the latter's interpola- 
tion operator simply sets SM in a child cell to be 
1 / r 3 that of the parent cell (for a three-dimensional 
spatial grid). 

(v) Following injection of arrays D and SM from level 
L+l to level L, where the injection operator for SM 
is an algebraic sum over child cells (a) zero all type 
B cells in SM on level L + l, (b) call the application 
routine (first item in the next list) to apply the mass 
corrections to D stored in the injected SM to type 
A cells on level L, (c) zero all type A cells in SM 
on level L. 



In other words, an inner (non-physical) boundary on level L must 
be at least one cell away from any inner boundary on level L — 1 . 
If the hierarchy is generated by truncation error which is suffi- 
ciently smooth, inner boundaries will typically not be coincident. 
Also, experience suggests it is often more challenging to get an 
AMR evolution stable if inner boundaries are too close, so in all 
this restriction is not particularly limiting. 



The following are new tasks that the unigrid application 
code needs to implement: 

(i) A routine that will add the mass corrections stored 
in SM to D for all type A cells on a given grid (i.e., 
set Dl — > Dl + SM/V L 13. 

(ii) When taking a single time step on a grid, for any 
cell marked type A, set SM to minus the change 
in mass of the cell from fluxes through cell faces 
indicated by the bitmask. For example, with the 
case illustrated in Fig |T] and discussed above around 
Eqs. (HU) and O, set SM L = -V L SD L . 

(iii) When taking a single time step on a grid, for 
any cell marked type B, add to SM the change 
in mass of the cell from fluxes through cell faces 
indicated by the bitmask. For example, with the 
same example above, set 8ML+x,j,k — > SM l + \ + 

VL+l,j,k, m SDL+l,j, 

For the GR-hydro equations we have five conserved 
fluid variables, D and S a . Though the latter do have 
nonzero source terms — since gravity can be a source 
(or sink) of energy-momentum — the above algorithm 
ensures there will be no artificial loss/gain in the pres- 
ence of AMR boundaries due to truncation error from 
the advection terms. 



III. TESTS 

In this section we present a number of tests of the 
methods presented above. We begin by demonstrating 
the fourth-order convergence of the evolution of the Ein- 
stein equations for vacuum spacctimes before moving on 
to a number of flat space, rclativistic hydrodynamics tests 
that probe the treatment of fluid discontinuities. Wc 
conclude with several tests of hydrodynamics in curved 
spacetimes. 

A. Vacuum evolution 

In [ill, [HI several tests of convergence of an earlier 
version of the code (without hydrodynamics) were pre- 
sented. However, since then we have updated the evo- 
lution of the Einstein equations to fourth-order spatial 



Since wc consider D a density and 6M a mass, this requires 
normalization by the volume clement Vl, which the application 
knows. Note that in our code even though we have included 
the uncompactified metric volume clement \J—g in the definition 
of the conservative variables and fluxes, compactification (and 
in axisymmetry, the cylindrical shell volume element) effectively 
makes the grid non-uniform and so the volume scaling is non- 
trivial. An alternative implementation could move this correc- 
tion step to the driver list of tasks, though then the application 
would need to supply the driver with the array of local volume 
elements. 
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differencing and fourth-order Runge-Kutta time differ- 
encing, so we first show two vacuum tests: a Briil wave 
evolution [92, 9:| and a boosted BH evolution. 



1. Brill wave 

For the Brill wave test we begin with initial data where 
the spatial line element is given by 



ds 2 = ^(e B dx 2 + eBy \ +z2 dy 2 
(e * ~ l)VZ (dydz + dzdy) + eBz2 + y2 dz 2 ) (22) 



where r = \Jy 2 + z 2 , B = 2Ar 2 exp(—(r/a r ) 2 — (x/a x ) 2 ), 
and the value of the conformal factor is determined by 
solving the Hamiltonian constraint. We choose A = 40, 
ay = 0.16, and a x = 0.12. The initial data is chosen to be 
time symmetric (7^ = 0) and maximally sliced (K = 0) 
with the conformal lapse a := ^>~ 6 a = 1. The remaining 
metric components are chosen to satisfy the harmonic 
gauge (Dx a = 0). This describes a gravitational wave 
that initially collapses inward before dispersing. In Fig. [2] 
we show results from convergence tests in axisymmetry at 
three resolutions where the medium and high runs had, 
respectively, 1.5 and 2 times the resolution of the low 
run. The constraint equations (H a — Ox a = 0) as well as 
the metric components show the expected fourth-order 
convergence. 



2. Boosted BH evolution 

As an additional vacuum spacetime test we evolved a 
boosted BH in three dimensions. We began with initial 
data describing a BH in harmonic coordinat es 19411 with 
boost parameter v — 0.25. As described in |42j, during 
the evolution we avoid the BH singularity by searching 
for an apparent horizon and excising a region within. 
To demonstrate convergence we performed this simula- 
tion at three resolutions, the lowest of which has approx- 
imately 30 points covering the diameter of the BH. The 
medium (high) resolution has 1.5 (2.0) times the number 
of points in each dimension, respectively. For all res- 
olutions we used the same AMR hierarchy, determined 
based on truncation error estimates at the lowest reso- 
lution, with six levels of 2:1 refinement. In Fig. [3] we 
demonstrate that the constraint equations are converging 
to zero at fourth-order. When hydrodynamics is included 
the theoretical limiting convergence rate of our code will 
drop to second-order (in the absence of shocks). However 
in vacuum dominated regions, for example the gravita- 
tional wave zone, one can expect that for the finite res- 
olutions we can practically achieve the convergence will 
be somewhere between second- and fourth-order. 
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FIG. 2. Top: The natural log of the L 2 norm of the con- 
straint violation, C a '■= H a — Ox a , for a Brill wave evolution 

natural log of <\J J \C a \ 2 d 2 x/ J d 2 x). The three reso- 



lutions shown are scaled assuming fourth-order convergence. 
Time is shown in units of M, the total ADM mass of the space- 
time, and the constraints are multiplied by M to make them 
dimensionless. The lowest resolution has a grid spacing of 
h — 1.56A/. Middle/Bottom: The value of the metric com- 
ponent gtt evaluated at (x,y,z) — (0, 50M, 0) (middle) and 
the difference in this quantity between low and medium reso- 
lution, and medium and high resolution (bottom), the latter 
scaled so that the two curves should coincide for fourth-order 
convergence. 



B. Relativistic hydrodynamic tests in flat 
spacetime 

We have performed a number of standard tests for rel- 
ativistic, inviscid hydrodynamics that probe how well a 
given numerical scheme handles the various discontinu- 
ities that arise. The best combination of reconstruction 
and flux calculation methods depends on the problem 



10 



x 10" 



m 



^ A I 

o 4 ; 







■ t ** »" ~ *:» - ~ * - V ^ 

> is -V >-* - • c > 






r 
» 

i 




i 
f 

i 






— Low res. 

— Med. res.x 1.5 4 
High res. x 2 4 



10 

time (M BH ) 



15 



20 



25 



FIG. 3. The L 2 -norm of the constraint violation (C a '■= 
H a — Di ) in the equatorial plane for a boosted BH simula- 
tion with v = 0.25. The three resolutions shown are scaled 
assuming fourth order convergence. Time is shown in units 
of A/bh, the ADM mass of the BH in its rest frame, and the 
norm of the constraints is multiplied by Mbh so as to make 
it dimensionless. 



under consideration. We have thus implemented several 
options and maintained a modular code infrastructure 
so that they are readily interchangeable and upgradable. 
While strong shocks such as the ones considered here are 
not expected to play an important dynamical role in bi- 
nary BH-NS mergers, they might be important in other 
potential applications of interest (such as NS-NS grazing 
impacts, or understanding EM emission from collisions). 
Thus, the ability to tailor the reconstruction and flux 
methods to the problem at hand may prove important 
in the future. In this section, we closely follow the se- 
quence of tests used in the development of the RAM code 
of Zhang and MacFadyen [95[, so that our results may 
be compared with theirs. Though they focus on more so- 
phisticated flux-reconstruction algorithms, their simpler 
methods (labeled U-PPM and U-PLM, denoting recon- 
struction of the unknowns with piccewisc parabolic and 
linear reconstruction, respectively) are comparable to the 
ones we employ. 



1. ID Riemann problems 

We first present a series of four relativistic, one- 
dimensional (ID) Riemann problem tests for which the 
exact solution is known (see Sections 4.1-4.4 of [Hj|). In 
all cases, the domain is a; € [0,1] and there are initially 
two fluid states, a left and a right, initially separated by 
an imaginary partition at x = 0.5. At t = 0, the parti- 
tion is removed and the fluid evolves to some new state. 



A T-law EOS is used for all the tests. In Table Q] we 
summarize the initial states and adiabatic indices used 
for the four tests, which we label as RT1 (Riemann Test 
1), RT2, RT3, and TVT (Transverse Velocity Test). We 
compare the performance of the various combinations of 
reconstruction schemes and flux methods to the exact so- 
lution and summarize the errors and convergence rates in 
Table [II] Exact solutions to these four tests were gener- 
ated using a solver provided by B. Giacomazzo, which is 
described in [96|. Taking HLL as our basic flux method, 
we performed this series of Riemann problem tests with 
four reconstruction methods: MC, minmod, WEN05, 
and PPM. For MC and WEN05, we also explored the 
effect of the flux method by running the tests with the 
Roe solver and Marquina's method. Most cases have a 
Courant-Fricdrichs-Lewy (CFL) factor of 0.5. However, 
the Roe solver, when combined with WEN05, does not 
seem to work well for problems with very strong shocks, 
such as RT2 and TVT. For a CFL factor of 0.5, we ob- 
tain acceptable results with Roe only by using a more 
diffusive limiter (like MC). For RT2 and TVT, we thus 
use Roe combined with WEN05 with a CFL factor of 
0.1. 

All of the methods we considered perform well on RT1, 
which is a fairly easy test. The lowest overall error occurs 
for WEN05 reconstruction (though the density profile 
between the shock and the contact discontinuity seems 
not to be as flat as in the other cases). The overall suc- 
cess of WEN05 may be due to the fact that the shock is 
relatively mild and there is an extended rarefaction that 
benefits from the high-order reconstruction. In Fig. 2] we 
compare the density profile obtained using HLL and var- 
ious reconstruction methods to the exact solution. We 
note that the tests which used the Roe or Marquina flux 
calculation with WEN05 do not have the oscillation vis- 
ible in the plot around x = 0.8 in the HLL-WEN05 case. 

The second Riemann test (RT2) is more difficult than 
the first, with the blast wave resulting in a very thin shell 
of material bounded by a shock on the right and a contact 
discontinuity on the left (see Fig. [5]) . The average con- 
vergence rates for this test show a marked difference be- 
tween the piecewise-linear and higher-order reconstruc- 
tion methods. WEN05 seems to perform best, but there 
is not much difference between HLL and Marquina or the 
Roe solver (with diminished CFL factor) with WEN05. 
As in RT1, the reconstruction method seems to be more 
important to the solution than the flux scheme. 

RT3 is a challenging problem in which the fluid on 
the left collides with the initially stationary fluid on the 
right, resulting in two shocks separated by a contact 
discontinuity. Our numerical solutions suffer from sig- 
nificant oscillations (particularly in the reverse shock) 
for all reconstruction schemes except PPM, which was 
specifically designed to suppress such post-shock oscilla- 
tions (see Fig. [6]). PPM also has the best convergence 
properties (0.85-1.16), with an average rate close to the 
expected value of unity. (Finite-volume hydrodynamic 
schemes such as this should converge at first order to a 
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Test T a P L 


PL 


Vl 


Pr 


PR 


VR 


RT1 5/3 13.33 


10.0 


0.0 


10" 8 


1.0 


0.0 


RT2 5/3 1000.0 


1.0 


0.0 


10" 2 


1.0 


0.0 


RT3 4/3 1.0 


1.0 


0.9 


10.0 


1.0 


0.0 


TVT 5/3 1000.0 


1.0 


0.0 


1.0 


1.0 (0.0,0.99) b 



TABLE I. The initial left and right states for the ID Riemann problems. 
a Adiabatic index of EOS 

b In this case v x =0 but the transverse velocity v y — 0.99 is nonzero. 



Reconstruction Flux method 




RT1 




RT2 




RT3 




TVT 






Error" 


Convergence 6 


Error Convergence 


Error Convergence 


Error Convergence 


MC 


HLL 


0.034 


0.82 


0.110 


0.59 


0.062 


0.77 


0.238 


0.72 




Roe 


0.032 


0.82 


0.110 


0.60 


0.052 


0.80 


0.233 


0.72 




Marquina 


0.036 


0.82 


0.127 


0.59 


0.056 


0.79 


0.227 


0.76 


Minmod 


HLL 


0.061 


0.86 


0.169 


0.42 


0.054 


0.71 


0.395 


0.76 


WEN05 


HLL 


0.033 


0.84 


0.093 


0.76 


0.039 


0.61 


0.191 


0.83 




Roe 


0.032 


0.85 


0.096 


0.79 


0.039 


0.60 


0.198 


0.81 




Marquina 


0.036 


0.85 


0.093 


0.76 


0.038 


0.66 


0.183 


0.82 


PPM 


HLL 


0.041 


0.88 


0.133 


0.67 


0.024 


1.01 


0.248 


0.78 



TABLE II. ID Riemann test results. 

"The LI norm of the error for resolution TV = 400. 

b The average convergence rate between runs with N = 200, 400, 800, and 1600. The ideal rate is unity for problems such as 
these containing discontinuities. 



weak solution of the equations when discontinuities are 
present.) 

For the transverse velocity test, the initial data are 
set up as in RT2, except that there is a transverse ve- 
locity v y = 0.99 on the right side of the partition. The 
strong shock propagates into the boosted fluid, and the 
structure of the shock is altered, since the velocities in 
all directions are coupled through the Lorentz factor [98{ . 
Again, the reconstruction technique influences the result 
more than the flux calculation. For WEN05 reconstruc- 
tion, the errors for HLL, Roe, and Marquina are all very 
close in magnitude. WEN05 and PPM yield the best 
results overall. In Fig. [7] we show the density profile at 
different resolutions for HLL combined with WEN05. 



2. ID Shock-heating problem 

We next consider a one-dimensional shock-heating 
problem as in 



951 ] . which tests a code's conservation of 



energy as well as the ability to handle strong shocks. 
For this problem, the computational domain is x € [0, 1] 
with a reflecting boundary at x = 1. The fluid moves 
toward this boundary with an ultrarelativistic initial ve- 
locity oiv = 1 — 10~ 10 . The fluid has an initial density of 
p = 1.0 and a very small amount of specific internal en- 
ergy, e = 0.003. The EOS is a gamma-law with V = 4/3. 



When the relativistic fluid slams into the wall, its kinetic 
energy is converted into internal energy behind a shock 
which propagates to the left. Because the fluid is initially 
cold, essentially all of the heat is generated through this 
conversion. As explained in [95| . the shock speed and 
the compression ratio of the shock (or equivalcntly, the 
post-shock density) is known analytically. We evaluate 
our errors by calculating the LI norm of the density er- 
rors on the entire computational domain. The average 
rate of convergence is also calculated using this measure 
of error. 



We performed this test using HLL with five different 
reconstruction methods at four different resolutions (200, 
400, 800, and 1600 zones). Results are shown in Ta- 
ble [Hi] We find that, due to the extremely strong shock, 
there is a tendency for post-shock oscillations to form 
with less diffusive reconstruction schemes (see Fig. [8]). 
The WEN05 solution is afflicted with severe post-shock 
oscillations and exhibits poor convergence to the exact 
compression ratio. Very diffusive reconstruction schemes 
(zero slope and minmod) are comparatively quite success- 
ful and converge rapidly to the exact compression ratio. 
PPM, with its flattening step, gives the best convergence 
rate overall. 
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FIG. 4. Density at t = 0.4 for Riemann Test 1 (RT1) with 
different reconstruction methods and the HLL flux scheme at 
resolution N = 400. The inset shows the shock and contact 
discontinuity. The exact solution was generated using the 
code of [H. 



FIG. 6. Density at t = 0.4 for Riemann Test 3 (RT3), a 
collision problem, for different reconstruction methods with 
HLL at resolution N = 400. The post-shock oscillations are 
largest in the MC case and smallest for PPM. 
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FIG. 5. Density at t = 0.4 for Riemann Test 2 (RT2) at dif- 
ferent resolutions for HLL-WEN05. The average convergence 
rate in this case is 0.76. The thin shell of material between 
the shock and the contact is particularly difficult to resolve. 



FIG. 7. Density at t = 0.4 for the transverse test (TVT) at 
different resolutions for HLL-WEN05. The average conver- 
gence rate in this case is 0.83. 



3. Emery step problem 

Next we consider the two-dimensional (2D) Emery step 
problem [sTl, |99| , with the setup as in [95| . In this sce- 
nario, a fluid flows through a wind tunnel at relativistic 
speed and hits a step, which is represented by a reflect- 
ing boundary condition. The computational domain is 
(x,y) £ [0,3] x [0,1] - [0.6,3] x [0,0.2] where the sub- 
tracted region represents the step. At the left boundary, 



inflow conditions are enforced (as in the initial data), 
while at the right, outflow conditions are enforced. All 
remaining boundaries are reflecting. The fluid is initial- 
ized with density p — 1.4, velocity v x — 0.999, and a 
r = 1.4 EOS. The pressure is set to 0.1534, giving a 
Newtonian Mach number of 3.0. 

Higher-order reconstruction methods seem to be essen- 
tial for this test problem. We find that the MC limitcr 
performs poorly, regardless of the flux method. Although 
the MC simulation is stable, the bow shock formed as 
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Reconstruction Error" Convergence 6 


ZERO 


950 


0.94 


Minmod 


801 


0.92 


MC 


1500 


0.85 


PPM 


824 


0.96 


WEN05 


1670 


0.53 



TABLE III. Shock -heating test results. For this test we also 
compare to zero slope reconstruction, labeled "ZERO." 
a The LI norm of the error for resolution iV = 400. 
b The average convergence rate between runs with N = 200, 
400, 800, and 1600. 



Shock Heating Test 
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tial fluid states in the lower/upper, left/right quadrants 
are 

(p,P,v x ,v y ) LL = (0.5,1,0,0) 
(p,P,v x ,v y ) LR = (0.1,1,0,0.99) 
(p,P,v x ,v y ) UL = (0.1,1,0.99,0) 
(p,P,v x ,v y f R = (0.1,0.01,0,0). 

In this simulation the lower-right and upper-left quad- 
rants converge on the upper-right quadrant creating a 
pair of curved shocks. We use a T = 5/3 EOS. In Fig. [TOl 
we show results from simulations using HLL or the Mar- 
quina flux method combined with WEN05 or the MC 
limiter. The first three panels are from runs with resolu- 
tion of 400 x 400 and a CFL factor of 0.5 and are com- 
parable to [95| and the references therein. Though the 
main shock features are captured by all of the methods 
we considered, oscillations arising from the curved shock 
fronts are evident in varying degrees. The fourth panel is 
similar to the first but contains a refined mesh in the cen- 
ter that has the same resolution as the other three panels, 
while the remainder of the domain has half the resolu- 
tion. Though the majority of the flow is thus effectively 
derefined, the principal features remain the same. This 
is despite the fact that the shocks must travel through 
or along refinement boundaries, and the numerical shock 
speeds differ slightly on either side of such boundaries 
due to the different truncation errors. 



FIG. 8. Density at t = 2.0 for the shock-heating test with 
different reconstruction methods and the HLL flux scheme. 
The inset focuses on the shock front. 



the fluid reflects off the step is distorted by large am- 
plitude post-shock oscillations. These propagate down- 
stream, rolling up the boundaries between the different 
solution regions. The higher resolution runs with MC 
also have these features, but at shorter wavelengths and 
lower amplitude. PPM and WEN05 reconstruction per- 
forms much better, and these results are shown at two 
resolutions i n Fi g. (This figure can be compared to 
those of HE [Tool.) The PPM results appear slightly bet- 
ter than WEN05 at a given resolution, likely because of 
the deliberate oscillation suppression in the PPM algo- 
rithm. 



4- 2D shock tube problem 

As an additional test of these algorithms' ability to 
propagate strong, multidimensional shocks we consider a 
2D shock tube test. The computational domain (x, y) € 
[0, 1] x [0, 1] is divided into four equal quadrants. The ini- 



C. Hydrodynamic tests in curved spacetime 

1. Bondi accretion 

As a first test of our code's ability to simulate rela- 
tivistic hydrodynamics in the strong field regime, we con- 
sider Bondi flow. We set up our initial conditions with 
a sta tiona ry solution to spherical accretion onto a black 
hole 
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We use Kerr-Schild coordinates for the black 
hole metric. In order to test our code's ability to converge 
to the correct solution we measured how the conserved 
density D differed from the exact solution as a function 
of time for three resolutions. The lowest resolution has 
a grid spacing of h — 0.078Mbh, while the medium and 
high resolutions have twice and 4 times the resolution re- 
spectively. As shown in Fig. [TTJ \\D — D cxact \\ converges 
to zero at second-order. For this test we tried both the 
MC and WEN05 limiters (with HLL for the flux calcula- 
tion) . Though both had similar levels of error and showed 
the expected convergence, WEN05 had larger errors at 
low and medium resolutions. This is probably because, 
at lower resolutions, the larger WEN05 stencil extends 
farther inside the black hole horizon where there is larger 
truncation error. 
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FIG. 9. Density contours (30 equally spaced in the logarithm) for the Emery step problem. The upper (lower) two plots 
show results for resolution 240 x 80 (480 x 160). For each resolution, the upper plot shows results for WEN05 reconstruction, 
and the lower for PPM. The respective minimum and maximum densities, (p m in,Pmax), are (1.0,1.0 x 10 2 ), (0.55,1.1 x 10 2 ), 
(0.82, 1.1 x 10 2 ), and (6.8 x 10~ 2 , 1.1 x 10 2 ). 



2. Boosted NS 



As an additional test of our evolution algorithm, we 
considered a single TOV star with a boost of v = 0.5, 
with astrophysically relevant EOS (the HB EOS of fill ]) 
and mass (1.35 M ). We performed a convergence study 
at three resolutions, the lowest of which has approxi- 
mately 50 points covering the diameter of the star. The 
medium (high) resolution has two (three) times the num- 
ber of points in each dimension, respectively. The AMR 
hierarchy is identical in all cases, with 7 levels of 2:1 
refinement, and was determined using truncation error 
estimates from the low resolution run. Figure shows 
that the constraint violations show the expected second- 
order convergence to zero. We also compared the perfor- 
mance of different reconstruction methods (though using 
the HLL flux method throughout). In Fig. Q21 we show 
the maximum density of the NS as a function of time 



for various reconstruction methods. Though the drifts 
and oscillations in density converge away for all meth- 
ods, we find that WEN05 gives the least density drift 
compared to MC and PPM at a given resolution. The 
drift in maximum density with PPM has to do with the 
way this particular implementation enforces monotonic- 
ity at ext rem a, which results in a loss of accuracy (see for 
example |l02j |). Modifying the way the algorithm handles 
smooth extrema can reduce this effect. We implemented 
one such modification (Eqs. 20-23 from jl02j ). the results 
of which are labeled 'PPM alt.' in Fig. [H 



3. Boosted NS flux correction test 

As a demonstration of the flux correction algorithm 
(outlined in Sec. Ill D|) to enforce conservation across 
AMR boundaries, we perform an additional boosted NS 
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FIG. 10. Density contours (30 equally spaced in the logarithm) for the 2D Riemann problem with, from left to right, top 
to bottom: HLL and MC, HLL and WEN05, Marquina and MC, and HLL and MC with mesh refinement. The respective 
minimum and maximum densities, (p min , pma X ), are (1.1 x 10~ 2 ,7.0), (8.2 x 10~ 3 ,8.1), (9.1 x 10~ 3 ,7.1), and (7.6 x 10~ 3 ,7.0). 
For the first three simulations a resolution of 400 x 400 was used. For the final simulation, a refinement region (red box) was 
placed in the middle with equivalent resolution, while the remaining grid has half the resolution (i.e. this simulation has lower 
resolution overall). A CFL factor of 0.5 was used throughout. 
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FIG. 11. The L 2 norm of the difference between the numer- 
ical and exact value of the fluid quantity D (divided by the 
norm of the exact solution) for Bondi accretion with MC and 
WEN05. The low resolution was run with a grid spacing of 
h — 0.078Mbh while the medium and high has twice and 4 
times the resolution respectively. The results are scaled for 
second-order convergence. 



test. We use the same conditions as the low resolution 
simulation outlined above in Sec. lIII C~2l but with a differ- 
ent AMR hierarchy. In particular, we keep the hierarchy 
fixed so that the boosted NS will move to areas of succes- 
sively lower refinement. In Fig.Q31 one sees that without 
flux corrections there is a « 0.1% loss in fluid rest-mass 
as the NS moves off the highest refinement level, and 
a pa 0.8% loss as the NS moves off the next to highest 
refinement level. This change in the conserved fluid rest- 
mass comes from the fact that there is a slight mismatch 
in fluxes at the mesh refinement boundaries due to trun- 
cation error. With the flux correction routine activated, 
this error is eliminated, and the only change in the total 
rest-mass is due to the density floor criterion (i.e., the 
numerical atmosphere). As an indication of how the use 
of flux corrections affects energy and momentum we can 
also compare the integrated matter energy and momen- 
tum as seen by a set of Eulerian observers. The matter 
energy density is given by T ab n a ni, and the momentum 
density is given by pi = —Tfn a where n a is the timclikc 
unit normal to the constant t slices. These quantities in- 
volve combinations of the conserved fluid variables and 
the metric and are subject to truncation error, especially 
since in this simulation the NS is allowed to move to 
lower resolution. In addition, these quantities can vary 
with time due to gauge effects (though in this case, the 
variation due to gauge effects is subdominant to the vari- 
ation mentioned below) , so we use them as an indication 
of the effect of flux corrections by comparing them for the 
simulations with and without flux corrections to a simi- 



lar simulation where the AMR tracks the NS and there is 
therefore essentially no flux across AMR boundaries. At 
the end of the simulation the run without flux corrections 
has 0.91% less energy compared to a similar simulation 
where the AMR tracks the NS, while the run with flux 
corrections has only 0.13% less energy. The run without 
flux corrections has 4.0% less momentum (in the boost di- 
rection) while with flux corrections the comparative loss 
is 2.5%. 

This test is somewhat artificial since we have deliber- 
ately prevented the AMR algorithm from tracking the 
NS. As long as the AMR algorithm can keep the bound- 
aries away from areas of non-negligible flux (as it does 
when following the boosted NS in the test described in 
Sec. IIII C 2p the effect of the flux correction algorithm is 
small, at the level of the numerical atmosphere that gets 
pulled along with the star. However, in astrophysical ap- 
plications, situations may generically arise in which fluid 
crosses AMR boundaries. For example, the tidal tails 
formed by the disruption of a NS by a BH will cross 
refinement boundaries, and likewise for the subsequent 
accretion disk that forms, since it would be much too 
costly to keep these entire structures on the finest mesh. 
Of course, the hydrodynamic solution is still subject to 
truncation error, which could in principle affect aspects 
of the dynamics at the same order of magnitude as pu- 
tative nonconservative effects. Though for certain prob- 
lems, such as calculating the amount of unbound mate- 
rial following a BH-NS merger, or studying the late time 
accretion, it could be quite advantageous to ensure con- 
servation within the hydrodynamic sector. It would be 
an interesting computational science problem to system- 
atically study the efficacy of AMR boundary flux conser- 
vation in such scenarios. 

Finally, we note that additional convergence test re- 
sults from this code were presented in [69} for the partic- 
ular BH-NS merger simulations discussed there. 



IV. CONCLUSIONS 

Numerous scenarios that fall within the purview 
of general relativistic hydrodynamics are still mostly 
unexplored — especially CO mergers involving neutron 
stars. There is a rich parameter space, of which large 
areas remain uncharted due to uncertainty or potential 
variability in BH and NS masses, BH spin and align- 
ment, the NS EOS, and other aspects. Beyond the pure 
hydrodynamics problem, the roles of magnetic fields and 
neutrino physics are just beginning to be explored by 
various groups, and we expect to add support for such 
physics to our code in the future. The potential appli- 
cations of robust and flexible numerical algorithms for 
evolving hydrodynamics together with the Einstein field 
equations are manifold. With this in mind, we have im- 
plemented methods for conservatively evolving arbitrary 
EOSs, in particular for converting from conserved to 
primitive variables without knowledge of derivatives; and 
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FIG. 12. The L 2 -norm of the constraint violation (C a := 
H a — Ox a) in the equatorial plane for a boosted NS simula- 
tion with v — 0.5 (using HLL flux calculation and WEN05 
limiter). The three resolutions shown are scaled assuming 
second-order convergence. Time is shown in units of Mo, the 
ADM mass of the NS in its rest frame, and the norm of the 
constraints is multiplied by Mq so as to make it dimensionless. 



FIG. 13. The relative variation of the maximum central 
density from its initial value (p m ax/pmax(t = 0) — 1) during 
a boosted NS simulation with v — 0.5 for various reconstruc- 
tion methods and for two different resolutions in the case of 
WEN05. 



we have implemented numerous reconstruction and flux 
calculation methods that can be used interchangeably to 
meet problem specific requirements. Though accurate 
treatment of shocks may not be crucial for BH-NS merg- 
ers (where shocks are not expected to be dynamically 
important), the same is not true of NS-NS binaries, espe- 
cially eccentric ones where the stars may com e into con- 
tact during nonmerger close encounters [l03| |. We have 
also taken care to implement a flux correction algorithm 
that preserves the conservative nature of hydrodynamical 
advection across AMR boundaries. Though strict conser- 
vation is not, strictly speaking, essential (since any non- 
conservation would be at the level of truncation error), it 
is an especially appealing property when studying, for ex- 
ample, CO mergers as potential SGRB progenitors. Af- 
ter merger, material that did not fall into the black hole 
— typically on the order of a few percent of the original 
NS mass — will fill a large volume making up an ac- 
cretion disk and potentially unbound material. Though 
accurately tracking this material is not important for the 
gravitational dynamics, it is critical for characterizing po- 
tential EM counterparts to the merger. 
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Appendix A: Spectral decomposition of the flux 
Jacobian 

Our conservative formulation of the hydrodynamical 
equations (|3I4[) can be written in vector notation as 5*q+ 
di(F l ) — S where q is a five dimensional vector of the 
conserved (in the absence of sources S) fluid variables 
q = (D,S t ,S x ,S y ,S z ) T and the flux F ?; = (Dv\(S t - 
y/—gP)v l , SjV 1 + Sjy/ — gP) T , where the index j in the 
flux is shorthand for the 3 components (x,y,z). Some 
flux calculation methods such as the Roe solver [7]} and 
the Marquina flux [78| require the spectral decomposition 
of the Jacobian which we give here. (See f 1041 ] for 
I 



the spectral decomposition for a similar formulation with 
slightly different conserved variables.) The eigenvalues 
are 

A± = aq(a ± 6) - (3 l (Al) 

and 

A 3 = aW - ft (A2) 

(with multiplicity 3), where a = (1 — c 2 )U l , b = 
c s V(l - C/ 2 )[ 7 "(l - IP<*) - aU% q = (1 - U 2 cl)~\ c s 
is the sound speed, and a, (3 l , 7 y are metric components 
as in (TTI))) . Here and throughout we use i € {x 7 y,z} to 
refer to the direction of the flux in the Jacobian with 
which we are concerned, In the following equations 
we use the index j as a shorthand for the three spatial 
components of the eigenvectors (that is, the components 
associated with S Xl S y , and S z ). The indices I and m 
are fixed by i and the indices n and p are fixed by j as 
indicated below. The index k is the only index that is 
summed over. A set of linearly independent right eigen- 
vectors is given by 



r± = (l, hW[U k p k - {a{f l - U l U l ) + Aft}/B], hW(Uj - 8}A/B)) , (A3) 



where A = [U l c 2 (l — U ) =F % and B = 7" — U l (a ± b)q, where k 



op 



and, 



r 3 = (k/{HW{k/ P -c 2 s )), U k f3 k - a,Uj) T , (A4) 



r 4 = (WU h 2hW 2 (U k f3 k -a)U l +hf3 h h( 1:ll +2W 2 U J Ui] 

(A5) 

where for r 4 , / = y,z,x for i = x,y,z respectively. The 
expression for rs can be obtained simply by replacing I 
with m, where m = z,x,y for i = x,y,z respectively, 
in the above expression for r 4 . H and W are as defined 
following (fj"2j) . 

We also give the corresponding left eigenvectors. 
Component-wise, for 1± = 



Tf 



l 3 ± = Tf 



{K - l){- 7 ^ + V T (W 2 ti - T lm )} + KW 2 V T q /a 

hinj mp - 7i P 7m„){l - KA l - (2K - l)V T U 1 } + (2K - l)V T iW 2 W 



(A6) 



where I = y,z,x and m = z,x,y for i = x, y, z respec- 
tively, and n = y,z,x and p = z,x,y for j = x,y,z re- 
spectively and T im = JUlmm-yimllm, £ = Ti m -^U l U\ 

K = (l-c^p/n)- 1 , A± = (a±b)q, V± = (tf*-A±)/( 7 «- 
U l A±), A 1 = (7" - tr i Z7 i )/(7 << - U' l A T ), and 

f- 1 = 2hWbq£(K - 1)( 7 " - WW) x 



[(7" - U l K+){ 1 11 - WA-)}- 1 . 

Furthermore, 

1 3 = ^-(n- c 2 sP ) (h, W/a, W(W - fi/aj) , (A7) 
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and the components of I4 and I5 are 
if = 



l\ = G, m (ahO _1 



(A8) 



where G ir „ = (j mm Ui - -fi m U m ) and for 1 4 , l = y,z,x and 
to = z,x,y for i = x,y,z respectively. The expression for 
I5 can be obtained from the above expression for I4 simply 
by interchanging I and m. 
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